Modeling multipolar gravitational-wave emission from small mass-ratio mergers 
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Using the effective-one-body (EOB) formalism and a time-domain Teukolsky code, we generate 
inspiral, merger, and ringdown waveforms in the small mass-ratio limit. We use EOB inspiral and 
plunge trajectories to build the Teukolsky equation source term, and compute full coalescence wave- 
forms for a range of black hole spins. By comparing EOB waveforms that were recently developed for 
comparable mass binary black holes to these Teukolsky waveforms, we improve the EOB model for 
the (2, 2), (2, 1), (3, 3), and (4, 4) modes. Our results can be used to quickly and accurately extract 
useful information about merger waves for binaries with spin, and should be useful for improving 
analytic models of such binaries. 
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I. INTRODUCTION 

Since the numerical relativity breakthrough of 2005 [1- 
3] , there have been tremendous advances both in the com- 
putation of gravitational radiation from binary black-hole 
systems, and in analytical modeling of this radiation us- 
ing approximate techniques. Despite rapid and ongoing 
advances, it remains a challenge for numerical relativity 
to quickly and accurately compute models that span large 
regions of parameter space. Extreme conditions such as 
large spins and small mass ratios are particularly chal- 
lenging, although there has been excellent recent progress 
on these issues [4-6]. 

The effective-one-body (EOB) formalism [7-17] makes 
it possible to analytically model the three main phases 
of binary black-hole evolution: inspiral, plunge- merger, 
ringdown. EOB has been used to model the dynam- 
ics and gravitational-wave emission from comparable- 
mass binaries [18-28], extreme mass-ratio inspiraling bi- 
naries [29-31] (neglecting conservative self-force effects), 
and small mass-ratio non-spinning binaries [32-36]. In 
order to study the transition from inspiral to plunge- 
merger and ringdown, Refs. [32, 33] suggested combining 
EOB with black hole perturbation theory. Concretely, 
they used EOB in order to compute the trajectory fol- 
lowed by an object spiraling and plunging into a much 
larger black hole, and then used that trajectory to de- 
scribe the source for the time-domain Regge-Wheeler- 
Zerilli (RWZ) equation [37, 38] describing metric pertur- 
bations to Schwarzschild black holes. They were then 
able to compute the full RWZ coalescence waveform and 
to compare with the EOB model. This was used in 
Ref. [34] to produce gravitational modes beyond the lead- 
ing (2, 2) mode and compute the recoil velocity. More re- 
cently, Refs. [35, 36] have used information from the RWZ 
modes to improve the modeling of the subleading EOB 
modes. A particularly beautiful feature of Ref. [36] is the 
use, for the first time, of hyperboloidal slicings in such 



an analysis. This effectively compactifies the computa- 
tional domain so that waveforms at future null infinity 
can be read out of the numerical calculation with great 
accuracy. 

References [39-41] developed a time-domain compu- 
tational framework based on the Teukolsky equation 
[42], which describes curvature perturbations of rotating 
(Kerr) black holes. The goal of these papers has been to 
understand gravitational waves produced by physically 
reasonable but otherwise arbitrary trajectories of small 
bodies bound to rotating black holes (such as slowly in- 
spiraling orbits, or trajectories that plunge into the hole's 
event horizon). This has been used to understand the 
small-mass ratio limit of merging black holes, studying 
for example the dependence of recoil velocity on black 
hole spin in this limit. This Teukolsky code has been op- 
timized to make effective use of modern many-core pro- 
cessor architectures, such as Graphics Processing Units 
(GPUs) [43]. 

In this paper, we combine EOB with the time-domain 
Teukolsky code developed in Refs. [39-41] to extend the 
ideas of Refs. [32, 33, 35, 36] in several directions. Our 
primary extension is, for the first time, producing full 
coalescence waveforms describing inspiral, merger, and 
ringdown for quasi-circular equatorial orbits in the Kerr 
spacetime. The energy flux we use in the EOB equations 
of motion comes from the factorized resummed wave- 
forms of Refs. [13, 15]. For the Schwarzschild limit, we 
model analytically three subleading modes [(2, 1), (3,3), 
and (4,4)] plus the dominant (2,2) mode, finding useful 
information about the plunge- merger, which we use to 
improve the comparable mass EOB model described in 
Ref. [28]. For more general spins, we calibrate the lead- 
ing EOB mode for spins a/M = -0.9, -0.5, 0.5, and 0.7. 
We also extract some information regarding subleading 
modes, and regarding the high prograde spin a/M = 0.9. 
These results for spinning binaries provide valuable input 
for improving the spinning EOB model of Refs. [14, 44], 
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as well as the spinning EOB waveforms of Refs. [23, 45]. 
This will in turn make it possible to develop models that 
can cover a much larger region of parameter space, in- 
cluding higher modes and extreme spins. 

Several other groups have also been using perturba- 
tion theory tools recently to improve our understanding 
for comparable mass and intermediate mass ratio bina- 
ries. For example, in Refs. [6, 46, 47] the authors di- 
rectly employ a moving-puncture trajectory (or a post- 
Newtonian-inspired fit to it) in the RWZ equation. They 
compare the resulting RWZ waveform with the results 
of full numerical relativity calculation for mass ratios 
1/15 and 1/10, finding good agreement. Another recent 
suggestion is the hybrid approach of Ref. [48], in which 
inspiral-plunge intermediate-mass black-hole waveforms 
are computed by evolving the EOB equations of mo- 
tion augmented by the perturbation-theory energy flux. 
An important issue in all attempts to model binary coa- 
lescence with perturbation theory is the computation of 
the so-called excitation coefficients, or more generally the 
question of which fundamental frequencies contribute to 
the radiation. In this context, perturbation-theory cal- 
culations are offering new insights [49-52]. 

The remainder of this paper is organized as follows. 
We begin in Sec. II by reviewing the EOB formalism 
for a test particle moving along quasi-circular, equato- 
rial orbits around a Kerr black hole. We then describe 
(Sec. Ill) the time-domain Teukolsky equation calcula- 
tion we use to compute the gravitational radiation emit- 
ted from a test particle that follows our EOB-generated 
trajectory. This section discusses in some detail numeri- 
cal errors which arise from finite-difference discretization, 
and from the extrapolation procedure by which we esti- 
mate our waves at future null infinity. Since we began 
this analysis, the Teukolsky code we use has been up- 
graded to use the hyperboloidal layer method [53]. This 
upgrade came too late to be used throughout our analy- 
sis, but has been used to spot check our estimates of this 
extrapolation error. 

Our results are presented in Sec. IV. We begin by com- 
paring the leading and three subleading Teukolsky modes 
with a/M = to the corresponding EOB modes cal- 
ibrated to non-spinning comparable mass binaries [28]. 
We then improve the non-spinning EOB model by in- 
cluding some features we find in our test-particle-limit 
calculation. We next calibrate the leading (2, 2) EOB 
waveform with our Teukolsky-equation results for spins 
a/M = -0.9, -0.5, 0.5, and 0.7. We conclude our re- 
sults by discussing the challenges of calibrating sublead- 
ing modes and of modeling extreme spin configurations, 
such as ones with a/M > 0.9. Section V summarizes 
our main conclusions, and outlines some plans for future 
work. A particularly important goal for the future will 
be to move beyond equatorial configurations, modeling 
the important case of binaries with misaligned spins and 
orbits. 

Throughout this paper, we use geometric units with 
G = c=l. 



II. DYNAMICS AND WAVEFORMS USING 
THE EFFECTIVE-ONE-BODY FORMALISM 

The Hamiltonian of a non-spinning test-particle of 
mass n orbiting a Kerr black hole of mass M and in- 
trinsic angular momentum (or spin) per unit mass a is 



H = f3 l pi + a yV 2 + l ij Pi Pj > 



(1) 



where the indices i,j label spatial directions, and the 
functions introduced here are given by 
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(2a) 
(2b) 
(2c) 



t is the time index, and is the Kerr metric. Working 
in Boyer-Lindquist coordinates (t, r, 9, </>) and restricting 
ourselves to the equatorial plane 8 = ir/2, the relevant 
metric components read 



9 



A 

A 

1 / 4a 2 M 2 



fJ 



t<f> _ 



A V A 

2a M 
' rA ' 



+ r 



(3a) 
(3b) 
(3c) 
(3d) 



where we have introduced the metric potentials 

A = r 2 - 2Mr - a 2 , (4) 
A=(r 2 +a 2 ) 2 -a 2 A. (5) 

We replace the radial momentum p r with p r * , the mo- 
mentum conjugate to the tortoise radial coordinate r*. 
The tortoise coordinate is related to the Boyer-Lindquist 
r by 



r 2 +a 2 
ar = : dr . 



(6) 



Since p r diverges at the horizon while p r » does not, 
this replacement improves the numerical stability of the 
Hamilton equations 



dr A dH 



dt r 2 + a 2 dp r 



(r,p r *,Pj,) ■ 



Tt =M ^ = W, {r,1Pr ' ,n) 



(7a) 
(7b) 



dt _ r 2 +a 2 Or ^ r 'P^'W + ^ p ^ ' 

d P4> _ nK - 

dt 



(7d) 
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Our trajectory is produced by integrating these equations 
using initial conditions that specify a circular orbit. We 
typically find in our evolutions a small residual eccentric- 
ity on the order of 3 x 10 -4 . 

In Eqs. (7a)-(7d), radiation- reaction effects are in- 
cluded following the EOB formalism. For the <f> com- 
ponent of the radiation-reaction force we use the non- 
Keplerian (nK) force 



nK x- 



1 dE 



(8) 



where «o = (Mfl) 1 / 3 , and dE/dt is the energy flux 
for quasi-circular orbits obtained by summing over 
gravitational-wave modes (l,m). We use 



dE 

lit 16tt 



- 8 t 



(9) 



1=2 m=- 



The non-Keplerian behavior of the radiation-reaction 
force is implicitly introduced through the definition of 
hi m . To describe the inspiral and plunge dynamics, we 
use the modes 



insp- plunge _ , F 



(10) 



The coefficients A^ m describe effects that go beyond the 
quasi-circular assumption and will be defined below [see 
Eq. (17)]. The factors hj m are the factorized resummed 
modes, and are given by [13] 



lF _ i.We) A( e ) rp JS tm I n \t 
h l m - h lra b L lrnZ \Ptm) 



(11) 



Here, e = tt(£ + rn) is the parity of the multipolar wave- 
form. The leading term in Eq. (11), h\^ e \ is the New- 
tonian contribution 



Mv 



S^M^Y*--" , (12) 



hi' ' = n 

tm n 

where TZ is distance from the source, Y (0,<j>) are the 
scalar spherical harmonics, and the functions n^f m and 
ci +e (v) are given in Eqs. (4a), (4b) and (5) of Ref. [54] 
with v = fx/M. For reasons that we will explain in 
Sec. IV A, we choose 



v rn 



1 (i+e-2) 



(£,m) ? (2,1), (4, 4), (13a) 
(*,m) = (2,l),(4,4). (13b) 
The quantities and tq introduced here are defined by 



v$ ee Mfl r n ee MQ, \(r/M) 3/2 + a/M 
The function in Eq. (11) is given by 



2/3 



H(r,p r *,P4,) , e = 0, 



L 



P<p vn 



e = 1. 



(14) 



(15) 



The factor T( m in Eq. (11) resums the leading order log- 
arithms of tail effects, and is given by 

T(£ + l- 2imMQ) 



T(£ + l) 



^■n rn M Q g 2i m M SI log(2 mQ,r ) 



(16) 

where r = 2M/^. [54] and T(z) ee / °° i 2 " 1 ^* dt is the 
complex gamma function. The factor e lSem in Eq. (11) 
is a phase correction due to subleading order logarithms; 
Si m is computed using Eqs. (27a)-(27i) of Ref. [54]. The 
factor (pim) m Eq. (11) collects the remaining post- 
Newtonian terms, and is computed using Eqs. (29a)-(29i) 
and Eqs. (Dla)-(Dlm) of Ref. [54]. 

Finally, the function Ni m entering Eq. (10) is given by 

•V,„, = I i-a^ m 7%2- + ahi 




exp 



where the quantities a^" 1 and b^ em are non-quasicircular 
(NQC) orbit coefficients. We will explain in detail how 
these coefficients are fixed in Sec. IV. 

We conclude this section by describing how we build 
the final merger-ringdown portion of the EOB waveform. 
For each mode (£, m) we have 



^mergcr- 



-RD 



JV-1 
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n=0 



(*) = y ] A-Umn e 



-icr fmrl (t-t m ™ tch ) 



(18) 



where n labels the overtone of the Kerr quasinormal mode 
(QNM), N is the number of overtones included in our 
model, and Ai mn are complex amplitudes to be deter- 
mined by a matching procedure described below. The 
complex frequencies <Ji rnn = 0Ji mn — i/T[ mn , where the 
quantities LOi mn > are the oscillation frequencies and 
T £mn > are the decay times, are known functions of 
the final black-hole mass and spin and can be found in 
Ref. [55]. In this paper, we model the ringdown modes 
as a linear combination of eight QNMs (i.e., N = 8). 

The complex amplitudes Ai mn in Eq. (18) are deter- 
mined by matching the merger-ringdown waveform (18) 
with the inspiral-plunge waveform (10). In order to do 
this, iV independent complex equations needs to be spec- 
ified throughout the comb of width At: 



match ' 



Details on 

the procedure are given in Ref. [28]. The full inspiral(- 
plunge)-merger-ringdown waveform is then given by 



h 



tm 



, insp-plunge afj-tm ±\ , l 



tm 



'6(t-t 



match/ ' 



(19) 
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In this analysis, we focus on waveforms emitted by a 
test-particle of mass ijl orbiting a Kerr black hole. Thus, 
we shall set to zero terms proportional to v = fj,/M 
in Eq. (11), excepting the leading v term in Eq. (12). 
Throughout this paper we restrict ourselves to the case 
v = 1CT 3 . 



III. THE TIME-DOMAIN TEUKOLSKY CODE 

A. Overview: The Teukolsky equation and its 
solution 

The evolution of scalar, vector, and tensor perturba- 
tions of a Kerr black hole is described by the Teukolsky 
master equation [42], which takes the following form in 
Boyer-Lindquist coordinates: 



(r 2 
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IT 



(20) 



The coordinates, the mass M, the spin parameter a, and 
the function A are as defined in the previous section. The 
number s is the "spin weight" of the field. When s = ±2, 
this equation describes radiative degrees of freedom for 
gravity. We focus on the case s = —2, for which <!/ = 
(r — iacos9) 4: ip4> where 1(14 is the Weyl curvature scalar 
that characterizes outgoing gravitational waves. 

To solve Eq. (20), we use an approach introduced by 
Krivan et al. [56]. First, we change from radial coordi- 
nate r to tortoise coordinate r* [Eq. (6)], and from axial 
coordinate <j> to <j), defined by 

d4> = d(t>+ ^dr . (21) 

These coordinates are much better suited to numerical 
evolutions, as detailed in Ref. [56]. Next, we exploit ax- 
isymmetry to expand \& in azimuthal modes: 



</) m (t,r, 9) 



(22) 



This reduces Eq. (20) to a set of decoupled (2+1)- 
dimensional hyperbolic partial differential equations 
(PDEs). We rewrite this system in first-order form by 
introducing a momentum-like field, 



n, 



(r 2 



d r * 



Following Ref. [56], we set <\> m and Ii m to zero on the 
inner and outer radial boundaries. Symmetries of the 
spheroidal harmonics are used to determine the angular 
boundary conditions: For m even, we have dg<p m = at 
9 = 0, it; for m odd, <f) m = at = 0, 7r. 

The right-hand side (RHS) of Eq. (20) is a source term 
constructed from the energy-momentum tensor describ- 
ing a point-like object moving in the Kerr spacetime. The 
expression for T is lengthy and not particularly illumi- 
nating. For this paper, it is suffices to point out that T 
is constructed from Dirac-delta functions in r and 9, as 
well as first and second derivatives of the delta function 
in these variables. These terms have coefficients that are 
complex functions of the black hole's parameters and the 
location of the point-like object. Details and discussion 
of how we model the deltas and their derivatives on a 
numerical grid are given in Ref. [39]. The delta func- 
tions are sourced at the location of the point-like object; 
the source T thus depends on the trajectory that this 
body follows in the Kerr spacetime. In this analysis, we 
use a trajectory constructed using the EOB formalism to 
specify the small body's location. 

One point worth emphasizing is that the source term is 
scaled by a factor of \ ji [see Eq. (2.39) of Ref. [39]]; i.e., 
the source is inversely weighted by the rate of change of 
coordinate time per unit proper time experienced by the 
orbiting object. This means that the source term "red- 
shifts away" as the object approaches the horizon. As 
a consequence, when describing a body that falls into a 
black hole, the Teukolsky equation (20) smoothly tran- 
sitions into its homogeneous form, connecting the gravi- 
tational radiation from the last few orbital cycles to the 
Kerr hole's quasinormal modes in a very natural way. 
The same behavior is seen in other analyses which model 
plunging trajectories using black hole perturbation the- 
ory (e.g., Refs. [34, 49, 57]). 

We implement this numerical scheme with a Fortran 
code, parallelized using a standard domain decomposi- 
tion (on the radial coordinate grid), and with OpenMPI 1 
enabled message passing. Good scaling has been ob- 
served for several hundred processor cores. In this anal- 
ysis, we used 128 processor cores for computing each m- 
mode for all cases we studied. 



B. Waveforms and multipole decomposition 

Far from the black hole, ip4 is directly related to 
and hx via 



^4 = 



1 (d 2 h A 



.d 2 h y 



ld 2 h 
2~W 



(24) 



2 V dt 2 dt 2 

The waveform h = h+ — ih x is then found by integrating 
(23) V'4 twice, choosing constants of integration so that h 



where S 2 = (r 2 



,2\2 



—a 2 A sin 2 9. We then integrate this 



system using a two-step, 2nd-order Lax-Wendroff finite- 
difference method. Details are presented in Refs. [39, 40]. 



1 The open source version of MPI, the Message Passing Interface: 
http:/ / openmpi.org. 
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at very late times (long after the system's waves have 
decayed to zero). 

As detailed in Sec. Ill A, our computation naturally 
decomposes the field \& (and hence ^4 and the waveform 
h) into axial modes with index m. For comparison with 
EOB waveforms, it is necessary to further decompose into 
modes of spin-weighted spherical harmonics. Following 
standard practice, we define 



h= — h im (t, r)_ 2 Y em (6, < 



(25) 
(26) 



In these equations, -^Yim is a spherical harmonic of spin- 
weight —2. Defining the inner product 



{Y tm \f)= / dfi- 2 17m(M)/ 



(27) 



(where * denotes complex conjugation), extracting Ci m 
and hf m is simple: 



Ctm(t,r) 
hi m (t,r) 



K{Y em \i>i) 
K{Y lm \h) . 



(28) 
(29) 



The complex wave mode hg m can also be obtained from 
Cg m by integrating twice, again choosing the constants 
of integration so that hi m — > at very late times. 



C. Numerical errors 

Our numerical solutions are contaminated by two dom- 
inant sources of error: Discretization error due to our 
finite-difference grid, and extraction error due to comput- 
ing and associated quantities at finite spatial location 
rather than at null infinity. 



Here, n — 2 is the ratio of grid spacing between the two 
resolutions. 

To estimate the remaining error in this extrapolated 
solution, we compare h 2 3 \ and h^\. Let us define 



Ah 



,(3) 
'2.5 



l (3) 
'1.5 i 



hfi - h 



2 - 5 1-1/n 3 



(31) 
(32) 



is a fourth-order estimate of the Teukolsky solution 
h, assuming that errors in h 



(3) 

2.5,1-5 



are third order. Defin- 



ing the amplitude \h\ and phase cj) as 



h = \h\ 



,itt> 



the amplitude error <5|/j|/|/i| and phase error 
S\h\ ^ f Ah\ 



are 



\h\ 



(33) 

(34) 
(35) 



Figure 1 shows discretization errors for several gravi- 
tational modes hg m extracted at r* — 950M. For this 
case, the large black hole is non spinning (a = 0). Am- 
plitude discretization errors are steady over almost the 
entire waveform, until very late times. In all cases, 
^l^l/l^l ^ a few x 1CP 3 . Similar behavior is observed 
for phase errors. For most modes, 5<j) < a few x 10~ 3 ra- 
dians over the coalescence. The highest (I, m) modes we 
consider approach 10 -2 radian error at the latest times. 
Because higher (£, m) modes require higher grid densities 
to be resolved, they tend to have larger discretization er- 



2. Extraction error 



1. Discretization error 

As discussed in Ref. [40], our time-domain Teukolsky 
solver is intrinsically second-order accurate. Since we 
compute our solutions on a two-dimensional grid in tor- 
toise radius r* and angle 9, we expect our raw numerical 
output to have errors of order (dr*) 2 , (dff) 2 , and (dr*d9). 
We mitigate this error with a variant of Richardson ex- 
trapolation, which we now describe. 

Consider waveforms generated at three resolutions: 
h[ 2) at (dr*,d9) = (0.064M, 0.2); h {2) at (0.032M, 0.1); 

and h$ at (0.016M, 0.05). Superscript "(i)" means the 
solution is zth-order accurate. We convert from second- 
order to third-order accuracy using [58] 



h {3) - h {2) 

"■1.5 — n \ 



(2) 



,(2) 



1-1/n 2 ' 



■,(3) 
L 2.5 



,(2) 



,(2) 



,(2) 



I /n 2 



(30) 



The code used in the bulk of this analysis extracts ^ 
(and derived quantities such as tp4 and h) at large but fi- 
nite radius. These quantities are more properly extracted 
at future null infinity. Although it has very recently be- 
come possible to extract waveforms at future null infinity 
(see Refs. [36, 53]), we did not have this capability when 
we began this analysis. Instead, following Ref. [59], we 
extract waveforms at multiple radii, and fit to a polyno- 
mial in 1/r. Again defining amplitude \h\ and phase <f> 
using Eq. (33), we put 

\h\(t - r*,r) = \h\ (0) (t -r*) + J2 >( 36 ) 

k=l 

<j>(t - r*,r) = ( o) (t -r*) + J2 ~ • ( 37 ) 

k=l 

The time t — r* is retarded time, taking into account the 
finite speed of propagation to tortoise radius r*; A is 
the order of the polynomial fit we choose. The functions 




|/i|(o)(t — r*) and 4>(o)(t — r*) are the asymptotic am- 
plitudes and phases describing the waves at future null 
infinity. 

We extract waveforms at radii r = 150M, 350A/, 
550M, 750M and 950Af. We then perform non-linear, 
least-squares fits for \h\nA and 4>if.) using the Levenberg- 
Marquardt method [60] to find the asymptotic waveform 
amplitudes and phases. Following Ref. [59], we use N = 3 
for the order of our fit, and estimate errors by comparing 
the fits for N = 3 and N = 2. 

Figure 2 shows the extrapolation errors we find for the 
same case shown in Fig. 1. For most of the evolution, ex- 
trapolation errors are smaller than discretization errors. 
In particular, the amplitude errors are at or below 10 -3 
for most of the coalescence; phase errors are at or be- 
low 10 -3 radians. Both phase and amplitude errors grow 
to roughly 10~ 2 very late in the evolution. Note that 
the largest errors in ip^ come at the latest times, when 
the waves have largely decayed away. In other words, the 
largest errors occur when the waves are weakest. Because 
we compute tp4 and then infer amplitude and phase, both 
amplitude and phase are affected in roughly equal mea- 
sure by these late time errors [see Eqs. (34) and (35)]. 
Our numerical errors appear to be of similar size to error 
estimates seen in related analyses (e.g., Ref. [36]). 

Finally, it is worth noting that, thanks to the hyper- 
boloidal layer method introduced to time-domain black 
hole perturbation theory in Refs. [36, 53], it will not be 
necessary to perform this extrapolation in future work. 
The codes will, to very good accuracy, compute the wave- 
form directly at future null infinity. Although this ad- 
vance did not come in time for the bulk of our present 
analysis, we have used it to check our error estimates 
in several cases. We find that our total numerical error 
estimates (discretization plus extrapolation error, com- 
bined in quadrature) is similar to the errors we compute 



using the hyperboloidal layer method 2 . This gives us 
confidence that our error estimates are reliable. 



D. Comparing time-domain and frequency-domain 
Teukolsky codes 

As a further check on the accuracy of our numerical 
Teukolsky-based waveforms, we compare time-domain 
(TD) waveforms computed using the techniques de- 
scribed here with frequency-domain (FD) waveforms 
[61, 62]. Since we only calibrate the higher-order modes 
in the EOB model for a = 0, we focus on that case here. 
We expect our conclusions to be similar for spinning cases 
since we use the same procedure to estimate errors in that 
case. 

As described in Sees. I and III A, in our analysis the 
source term for the TD waveforms [cf. Eq. (20)] de- 
pends on the EOB inspiral and plunge trajectory. For 
FD waveforms, by contrast, the source is built from a 
purely geodesic trajectory. This is because the FD code 
uses the existence of discrete orbital frequencies. For this 
analysis, we specialize further to circular-orbit equatorial 
geodesies, but allow these geodesies to evolve adiabati- 
cally using FD Teukolsky fluxes, as described in Ref. [63]. 
Previous work has shown that a self-consistent adiabatic 
evolution implemented with our FD code is in excellent 



2 It is worth emphasizing that codes which use the hyperboloidal 
layer method are much faster than those which use the extrapo- 
lation described here; we find a speedup of roughly ten (for the 
scale of the evolutions performed in the context of this work). 
Although it is gratifying that these extrapolations reliably im- 
prove our numerical accuracy, the substantial speed-up means 
that upgrading our method is worthwhile for future work. 
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agreement with the EOB model during the inspiral [29], 
and so it makes sense to compare TD and FD waveforms 
during this phase of the coalescence. It is also worth 
noting that FD waveforms can generally be computed 
to near machine accuracy using spectral techniques [64]. 
The only limitation on their accuracy is truncation of the 
(formally infinite) sums over multipoles and frequency 
harmonics. We can thus safely assume that the difference 
between TD and FD waveforms is only due to errors in 
the TD waves. 

To perform this comparison, we align the I — m — 2 
TD and FD waveforms by introducing time and phase 
shifts At and A<p that minimize the gravitational phase 
difference at low frequencies. More specifically, we choose 
At and A(f> in order to minimize 

f 2 [0™ (t) - 0™ (t + At) + A<f\ 2 dt , (38) 
Jti 



(£,m) 


«~ TD 


cjTD 
°4>lm 


5 \ h \i m 


8\h\Jm 




\h\tm 


\h\tm 


(2,2) 


7.71x10" 


-4 


1.27x10 


-3 


1.22x10 


-4 


8.20x10" 


4 


(3,3) 


1.18x10" 


-3 


2.22x10' 


-3 


1.95x10 


-4 


1.64x10" 


-3 


(2,1) 


4.05x10" 


-4 


1.28x10 


-3 


2.55x10 


-4 


1.03x10" 


-3 


(4,4) 


1.58x10" 


-3 


2.94x10' 


-3 


2.80x10 


-4 


2.80x10" 


-3 


(3,2) 


7.71x10" 


-4 


3.92x10" 


-4 


4.09x10" 


-4 


3.13x10" 


-3 



TABLE I. The phase difference and fractional amplitude dif- 
ference for various modes, averaged over the time interval 
ti — t\ [see Eq. (38)]. We compare the amplitude and phase 
error found by comparing TD and FD waveforms (columns 2 
and 4) with the TD errors we estimate using the techniques 
discussed in Sec. Ill C. In all cases but one [phase error for 
the (3, 2) mode] , our numerical error estimates are larger than 
those we find comparing the two calculations; in that single 
discrepant case, the errors themselves are particularly small. 
This is further evidence that our numerical error estimates 
are reliable. 



where t\ and ti are separated by lOOOAf and correspond 
to Mw 22 « 0.108 and a Mw 22 « 0.111, respectively. This 
low-frequency alignment is necessary for three reasons. 
First, the time coordinate of the TD waveform includes 
the effect of the extraction radii of the data used for 
the extrapolation; the FD waveforms are truly extracted 
at future null infinity. Second, the initial phases of the 
TD and FD trajectories are not necessarily the same, 
which introduces a phase offset between the two models. 
Third and last, as discussed in detail in Ref. [39], TD 
waveforms include an initial burst of "junk" radiation, 
which must be discarded. During that burst, the TD and 
FD trajectories may accumulate a small phase difference. 
We have found that small changes to t\ and t% do not 
significantly affect the alignment. 

Once At and A<f> are fixed, we have no freedom to intro- 
duce further time or phase shifts for the other modes. For 
instance, the difference between the FD and TD phases 



for the mode (£, m) is 



/FD— TD 



The fractional amplitude difference is 



ritiFD-TD 
5 \ h \hn 



\h\J°(t) 



\h\J°(t + At) 



(39) 



(40) 



The At and A</> used here are the ones that minimize 
(38). 

Table I compares <5</>™~ TD and S\h\^~ TI) /\h\ em with 
the errors computed using the techniques described in 
Sec. IIIC. In particular, we examine the averages of 
^Jm - ™ an d <5Njm _TD /N^"i over the alignment in- 
terval (ti,ta)j an d compare them to the averages over 
the same interval of the TD numerical errors discussed 
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(^pcak 


O \ / i r 

- ipcak)/^ 


a/M 


/ j_ 22 
l^pcak 


Q 

''peak 


(2,2) 




-2.99 


-0.9 




1.60 


(3,3) 




0.52 


-0.5 




-0.08 


(4,4) 




2.26 







-2.99 


(2,1) 




8.82 


0.5 




-7.22 


(3,2) 




6.25 


0.7 




-12.77 








0.9 




-39.09 



Mfi max 

0.090 
0.106 
0.136 
0.193 
0.236 
0.320 



TABLE II. The time difference (t^ k - tp eak ) / M between the 
peak of the Teukolsky modes' amplitude and of the orbital 
frequency, for a/M = 0. 



in the previous section. For this comparison, we average 
the sum (in quadrature) of discretization and extrapola- 
tion errors. We see that the difference between TD and 
FD is always within the TD numerical errors, except for 
the £ — 3, m = 2 mode. This mode is among the weak- 
est of those that we show in Table I, which makes the 
extraction procedure described in the previous section 
considerably more difficult. 



E. Characteristics of time-domain Teukolsky 
merger waveforms 

We now turn to the waveforms produced by the TD 
Teukolsky analysis and their general characteristics. Fig- 
ure 3 examines the behavior of the dominant TD modes 
[(£,m) = (2, 2), (3, 3), (4, 4), (2,1), (3, 2)] during plunge, 
merger and ringdown. We also show the orbital frequency 
of the EOB trajectory used to produce the TD data. 

For the non-spinning case (left panel), the peak of the 
(2, 2) amplitude comes slightly earlier than the peak in 
orbital frequency, while higher-order modes peak later. 



A summary of the time difference ip" ak 



peak 



between 



the peak of the Teukolsky mode amplitude and that of 
the EOB orbital frequency is shown in Table II. This 
difference can be as large as 6.25M for the (3, 2) mode, 
and even 8.82M for the (2, 1) mode. 

The situation is different in the spinning case. A trend 
we see is that as the spin a grows positive, higher-order 
modes become progressively more important. In the 
right-hand panel of Fig. 3, we show the amplitude of 
the eight strongest modes for a/M — 0.9. As noticed in 
Ref. [54] , modes with £ = m tend to increase more during 
the plunge. For example, the (5,5) mode is smaller than 
the (3, 2) mode during inspiral, but becomes larger dur- 
ing the plunge. The (6, 6) mode also grows quickly during 
the plunge (which roughly ends at the time of the peak 
of the orbital frequency) . This behavior is not surprising, 
in fact should become more and more pronounced as the 
spin a/M — > 1; modes with large multipole moments be- 
come as important as low-£ modes in that limit, at least 
for quasi-circular orbits [65, 66]. 



Table III shows the time difference (t 22 



peak 



ipeak)/^ 



between the peak of the (2,2) amplitude and of the or- 
bital frequency for the values of spin that we consider 
in this paper. As the spin grows, the orbital frequency 



TABLE III. The time difference (i^ - r* caS5 )/M between 
the peak of the Teukolsky (2, 2) amplitude of the orbital fre- 
quency, for various values of the spin a/M. Also shown is the 
value of the orbital frequency at that peak, Mf2 max - 



peaks later and later relative to the peak of the (2, 2) 
amplitude. This has important implications for model- 
ing the EOB merger-ringdown waveform, as we discuss 
in detail later in the paper. 

Another interesting feature of the Teukolsky wave- 
forms that we find is shown in Fig. 4. This figure 
shows the gravitational-wave frequency (denned as the 
time derivative of the phase) for the (2, 2) mode during 
plunge, merger, and ringdown, for several spin values. 
Notice the strong oscillations seen at late times (during 
the final ringdown) for spins a < 0. These oscillations 
grow as the spin decreases, and become very large for 
a/M = —0.9. We have verified that these oscillations 
(even in the a/M = —0.9 case) are not numerical arti- 
facts. We have found that they are insensitive to numeri- 
cal resolution and floating-point precision, that they also 
appear in the context of other plunging retrograde tra- 
jectories. We find that they are due to a superposition of 
the dominant (2, 2) QNM with the (2, -2) QNM, which 
is also excited during the plunge [34, 56]. 

In order to extract the relative amplitudes of these 
modes, we fit the gravitational-wave frequency in a region 
where the higher overtones of the (2, ±2) modes have 
decayed away and the waveform can be described by 



h(t) = ho(e' 



J°2-20 t+P 



(41) 



where C2±20 = ^>2±20 — «/ r 2±20 are the complex QNM fre- 
quencies [with overtone n — 0, as introduced in Eq. (18)]. 
The complex parameter h and real parameters a, (3 are 
left unspecified. From Eq. (41) one can then calculate 
the frequency as 5J[— ih(t)/h(t)]. Because h cancels in 
this expression, we are left with the parameters a and /3, 
which can be determined by numerical fitting [34]. We 
find that the relative excitation a of the (2, —2) modes 
goes from a ~ 0.005 for a/M — 0, to a w 0.03 for 
a/M = -0.5 and a w 0.46 for a/M = -0.9. This is 
nicely in accord with the growing strength of the oscilla- 
tions as a/M —> — 1 that is seen in Fig. 4. 

A possible reason why the (2, —2) QNM is strongly 
excited for large negative spins can be understood by ex- 
amining the particle's trajectory. When a < 0, the spin 
angular momentum is oppositely directed from the or- 
bital angular momentum. During the inspiral, when the 
orbit is very wide, the orbit's angular velocity is opposite 




FIG. 3. Amplitude of the dominant modes during plunge, merger and ringdown for a = (left panel) and a/M = 0.9 (right 
panel). We also show orbital frequency, scaled to fit on the plot. As expected, the orbital frequency asymptotes to the horizon's 
angular velocity at late times, because the frame dragging locks the particle's motion to that of the horizon. The vertical 
dashed line in the two panels marks the position of the peak of the orbital frequency. 




0.8 1 i t j 0.8 




° 5150 5200 5250 17050 17100 17150 
t/M t/M 

FIG. 4. The evolution of wave frequency for the (2, 2) mode 
for several spins. Notice that the strong late oscillations (dur- 
ing ringdown), which are particularly prominent as the spin 
decreases towards negative values; they are especially large 
for a/M = —0.9. These oscillations are physical, and arise 
from superposition of the (2, ±2) QNMs. 



to the sense in which the horizon rotates. At late times 
(during the final plunge), the particle's motion becomes 
locked to the horizon by frame dragging. The particle's 
angular velocity thus flips sign at some point during the 
plunge when a < 0. This change in angular velocity is 
most pronounced for large negative spins, since the dif- 
ference between the frequency at the innermost stable 
circular orbit (ISCO) and at the event horizon is largest 



for large negative a. 

Figure 5 shows, as an example, the EOB trajectory 
we used to produce the a/M = —0.9 Teukolsky wave- 
forms. As viewed here, the horizon rotates in the clock- 
wise sense. After the anti-clockwise inspiral, the parti- 
cle plunges and its angular velocity flips sign before the 
particles settles on a quasi-circular orbit with r — > r + 
and f2 — > fl+ = a/2Mr + as t — s- +oo (where r + = 
M + y/M? — a 2 is the coordinate radius of the event 
horizon). This behavior leads us to conjecture that the 
(2, —2) QNM is excited by the last part of the plunge, 
when the particle is corotating with the black hole. The 
(2, 2) QNM is excited by the final inspiral and initial 
plunge, when the particle is counter-rotating relative to 
the black hole. When a > the particle's motion is al- 
ways co-rotating with the black hole, both during inspi- 
ral and through the plunge. This conjecture thus ex- 
plains why oscillations in the ringdown frequency are 
much less significant for a > and seem to disappear 
when a/M s» 1 (see Fig. 4). 



IV. COMPARISON OF THE EOB MODEL 
WITH THE TEUKOLSKY TIME-DOMAIN 
WAVEFORMS 

In this Section, we present the main results of this 
paper, comparing the EOB waveforms for binary coales- 
cence with waveforms calculated using the time-domain 
Teukolsky equation tools described in the previous sec- 
tion. We begin by comparing Teukolsky waveforms (for 
a = 0) with an EOB model that has been calibrated for 
the comparable-mass case (Sec. IV A). The agreement is 
good for some modes [(2,2) and (3,3)], but is much less 
good for others [(2, 1) and (4, 4)]. We nail down the rea- 
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r cos(()))/M 



FIG. 5. The EOB equatorial orbit used to produce the 
Teukolsky waveforms for a/M = —0.9, focusing on the tran- 
sition between the inspiral and the plunge. As shown here, 
the Kerr black hole rotates clockwise. The red external circle 
is the ISCO; the blue internal one is the event horizon. The 
particle initially moves anti-clockwise on quasi-circular orbits, 
but flips to clockwise motion during the plunge as its angular 
motion becomes locked to the horizon's motion. 



son for this disagreement, recalibrate the EOB model, 
and show much better agreement in the comparison in 
Sec. IV B. We then consider a =^ 0. Focusing on the 
(2, 2) mode, we compare Teukolsky and EOB waveforms 
for a range of spins in Sec. IV C. 



We stress that all the comparisons that we present 
in this paper have been performed between EOB and 
Teukolsky waveforms produced with the same EOB tra- 
jectory. For example, in order to recalibrate the EOB 
model, we start with a reasonable EOB trajectory, we 
feed that to the Teukolsky code, and compare the result- 
ing Teukolsky waveforms to the EOB waveforms. If the 
waveforms do not agree, we modify the EOB trajectory so 
that the EOB waveforms agree with the Teukolsky wave- 
forms produced with the old EOB trajectory, and then 
feed the new EOB trajectory to the Teukolsky code. We 
then iterate until this procedure has converged. 




0.05 0.1 0.15 0.2 0.25 
V 



0.05 0.1 0.15 0.2 0.25 
V 



FIG. 6. The amplitude h (top panels) and gravitational-wave 
frequency oj (bottom panels) when the (2,2) (left) and (4,4) 
(right) modes reach their peak. Circles at v > 0.12 denote 
data points extracted from the numerical-relativity simula- 
tions; the left-most points at v — 10 -3 are data extracted 
with the Teukolsky code. The solid lines are quadratic fits to 
the data points. 



A. Comparison of comparable-mass EOB 
waveforms and Teukolsky waveforms for a — 

Reference [28] presents an EOB model calibrated to 
numerical-relativity simulations of non-spinning black- 
hole binaries with mass ratios mz/mi = 1, 1/2, 1/3, 1/4 
and 1/6. This model achieves very good agreement 
between the phase and amplitude of the EOB and 
numerical-relativity waveforms; see Sees. II and III of 
Ref. [28] for details. As background for the comparison 
we will make to the Teukolsky waveform, we briefly dis- 
cuss how the EOB inspiral-plunge waveform was built, 
and how the merger-ringdown waveform was attached to 
build the full waveform 



For each mode, Ref. [28] set a 4 
in Eq. (17), and fixed the remaining coefficients a. 
(1,2,3)] and 
five conditions: 



b h tm = 





i e 



£ (1,2)] by imposing the following 



1. The time at which the EOB /122 reaches its peak 
should coincide with the time at which the EOB 
orbital frequency f2 reaches its peak. We denote 
this time with £p Cak - The peaks of higher-order nu- 
merical modes differ from the peak of the numerical 
ft-22; we define this time difference as 



At 



Ira 
peak 



= t 



- 1 



22 



peak peak 
= td\tif£\/<tt=0 ' 



t. 



d\h^\/dt=0 ■ 



(42) 



We require that the peaks of the EOB h( m occur 
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FIG. 7. Comparison between (2,2) and (3,3) modes generated by the Teukolsky code for a = 0, and the corresponding modes 
produced with the EOB model of Ref. [28] . The EOB model was calibrated to numerical-relativity simulations at mass ratios 
1, 1/2, 1/3, 1/4, 1/6, and extrapolated to v = 10 -3 . Upper panels show the real part of the modes; lower panels show the phase 
and fractional amplitude differences. 



at the time t£ eak 



d\hf™\ 



dt 



^-"pcak- 



,+At£ 



(43) 



2. The peak of the EOB ht m should have the same 
amplitude as the peak of the numerical he m : 



'EOB fj. n 



peak 



I t NR/Jm \\ 
\ n lm V c poakJ | 



(44) 



3. The peak of the EOB he m should have the same 
second time derivative as the peak of the numerical 

him- 



d 2 \hf^\ 



dt 2 



d 2 \hZ\ 



dt 2 



(45) 



4. The frequency of the numerical and EOB h? m wave- 
forms should coincide at their peaks: 



£m Upcak + £it pcak) ~ W «m UpcakJ ■ l 4b ) 



5. The time derivative of the frequency of the numer- 
ical and EOB hg m waveforms should coincide at 
their peaks: 



^Fm B ( i pcak + A£p" ak ) — w]^(t p ™ k ) • (47) 



[Note that the quantities hf® B referenced in the above 
equations are the same as the quantities h™^ plungc de- 
fined in Eq. (10).] 



The functions Atz 



\hZ(C^)l d 2 \hZ\/dt 



,NR (+tm 



(t^ ak ), and <«(t £ ™ ak ) described in Ref. [28] were 



extracted from numerical-relativity and Teukolsky data, 
and approximated by smooth functions of the symmet- 
ric mass ratio v. Least-square fits for these quantities 
were given in Tabic III of Ref. [28]. These fits included 
information about the v = 10 -3 case from the analysis 
of this paper (which was in preparation as Ref. [28] was 
completed) . 

Since |^m(ip™ ak )| and d 2 \hi m \/dt 2 \ lm approach zero 

F peak 

in the test-particle limit, their input values at v = 10 
do not affect the least-square fits very much. For the 
(2,2) and (3,3) modes, the data points are very regu- 
lar. This is illustrated for the (2, 2) case in Fig. 6. The 
residues of the fit at v = 10~ 3 are very small, and the v- 
fits agree well with the values from the Teukolsky code at 
v = 10~ 3 . Unfortunately, this is not the case for the (2, 1) 
and (4, 4) modes. Figure 6 shows this for the (4, 4) mode. 
The |^44(ipoak)l data points do not lie on a smooth curve, 
and so the fit is intrinsically unstable. By minimizing 
the relative residue instead of the absolute residue in the 
least-square fit, we increase the weight on the v = 10~ 3 
data point and get a much better fit at low mass ratio, 
but at the cost of a much poorer fit in the comparable- 
mass regime. The situation is even worse for u>44(t / ^, ak ) , 
for which the data points for comparable masses have 
a rather irregular trend. These results emphasize the 
need for more accurate numerical-relativity data describ- 
ing the higher-order modes in order to smoothly connect 
these quantities from the test-particle limit to the equal- 
mass case. 



Once the coefficients a^ m and 6^ fm are known, we 
calculate /j^P-P lun se us j n g jt^ (10), and attach the 
QNMs using Eq. (19). We assume the following comb 



12 



6x10 
3xl0' 4 


-3xl0' 4 
-6xl0 4 
0.01 
0.005 


-0.005 
-0.01 



EOB Re(h 21 ) R/M 
Teuk-code Re(h ) R/M 



4600 



4800 



5000 



f ^\ J fc 



5<j) h (rad) 
5lhl/lhl 



4600 



4800 
t/M 



5000 




5150 5200 5250 



5150 5200 5250 
t/M 



2x10' 



-2x10 



EOB Re(h 44 ) R/M 
Teuk-code Re(h 44 ) R/M 






4600 


4800 


5000 


0.02 












it 


0.01 


— 


s 


.> 





- &|> h (rad) 
8lhl/lhl 




-0.01 








4600 


4800 


5000 




-2x10 



: /"' 




1 




y 












\ 




\ 




: \ 









0.0 



-0.5 



■1.0 



t/M 



t/M 



FIG. 8. The same as Fig. 7, but for the (2,1) and (4,4) modes. 



widths [28], 



A^match — 5 M . 



a .44 
^match 

and choose i^ tch 



= 9M. 



A^ 3 atch = 12M, 

A *match = 8 M , 



(48a) 
(48b) 



peak 



A^™ ak with A^ ak given 



in Eq. (42). It was found in Ref. [28] that, after cali- 
brating the EOB adjustable parameters and aligning the 
EOB and numerical waveforms at low frequency, the dif- 
ference between t^ cak and £p Cak is typically ~ 1M. Thus, 
Ref. [28] assumed that tp eak = t^ 2 cak and consequently 
set A£pp ak = 0. Following these findings, we attach the 
QNMs at t£ eak for the (2, 2) mode, and at t% eak + A* p ™ k 
for all the other modes. 

Figures 7 and 8 compare the leading modes generated 
by this EOB model with the modes generated by the 
time-domain Teukolsky code. We adopt the waveform 
alignment procedure used in Refs. [21, 22, 28, 54] and 
Sec. HID, aligning the waveforms at low frequency by 
minimizing 



[Mt)-<h(t- At )- A </>] 2 dt, 



(49) 



over a time shift At and a phase shift A</>. Here, 4>\{t) 
and 4>2(t) are the gravitational phases of the EOB and 
Teukolsky h 2 2- We chose t% — t\ = 1000M, and center 
these times when the orbital frequencies are low. We 
have verified that our results are insensitive to the precise 
location of this integration interval, provided that it is 
chosen during the inspiral phase. 

As expected from the discussion above, Figure 7 
shows that there is quite good agreement between the 
EOB and Teukolsky models for the (2,2) and (3,3) 
modes. In particular, the difference in both the am- 
plitude and the phase is quite small until the inspiral 
reaches the ISCO. This excellent agreement is due to 
the resummed-factorized energy flux [13] employed in the 



EOB equations of motion and waveforms (see previous 
studies [15, 32-36, 67, 68]). The amplitude disagreement 
during merger and ringdown is due to our procedure of 
attaching QNMs to the EOB waveform [see Fig. 3 and 
discussion around in Ref. [28]]. The accumulation of some 
phase difference during plunge will be discussed at the 
end of this section. 

Figure 8 shows that the agreement between the EOB 
and Teukolsky (2,1) and (4,4) modes remains excellent 
during the long inspiral, but is not very satisfactory dur- 
ing the merger and ringdown. For the (4, 4) mode, the 
EOB amplitude becomes too large toward merger. This 
is a consequence of the excessively large residue of the 
z/-nt for |/i 4 4(tpo ak )| at v = 1(T 3 . For the (2,1) mode, 
the EOB model of Ref. [28] fails to reproduce a reason- 
able merger waveform. This problem is related to the 
fact that the value of At 21 cak at v = 10~ 3 used in Ref. 
[28] to determine the i/-fit is too large (see Table II). The 
problem is deeper than this, however. In particular, the 
value is uncertain due to the (unusual) broadness of the 
Teukolsky (2,1) mode's peak (see Fig. 3). We shall see 
in the next section that to improve the agreement of the 
(2, 1) mode, we need a smaller value for AipJ, ak . 

An additional source of error arises from the procedure 
that was used to compute the NQC coefficients a^ 22 , 22 , 
and a h z 22 used in N 22 [see Eq. (17)]. In Ref. [28], these co- 
efficients were calculated by an iterative procedure using 
the five conditions discussed at the beginning of this sec- 
tion. These coefficients have small but non- negligible ef- 
fects on the EOB dynamics: through the amplitude I/122I, 
they enter the energy flux [see Eq. (9)] and thereby influ- 
ence the rate at which the small body spirals in 3 . This 



3 Note that the NQC coefficients a^ Cm of higher-order modes con- 
tribute much less to the energy flux and can be safely ignored in 
the dynamics of comparable- mass binaries [28]. 
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FIG. 9. The same as Fig. 7, but for the Teukolsky-calibrated EOB model described in Sec. IV B. 



iterative procedure increases by a factor of a few the com- 
putational cost of generating \i22- To mitigate this cost 
increase, Ref. [28] suggested replacing the iterative pro- 
cedure with ^-fits for a^ 22 , 22 and Og 22 . These fits were 
obtained using data for mass ratios 1, 1/2, 1/3, 1/4 and 
1/6. The EOB waveforms shown in Figs. 7, 8 are then 
generated using these ^-fits extrapolated to v = 10~ 3 . 
When comparing the fit and the true values of 

a'i 22 . a 2 22 

and a 3 22 at v = 1(T 3 , we find a non-negligible difference 
which is responsible for ~ 0.4 rad difference between the 
EOB and Teukolsky (2, 2) modes, and for ~ 0.6 rad dif- 
ference for the (3, 3) modes. In the next section, we shall 
show that by returning to the iterative procedure, rather 
than using the fits, we can do much better. 

Lastly, we comment on why for the (2,1) and (4,4) 

/rn- 



modes in Eq. (13b) we replaced v^ + ^ with v. 



,C«+e-2) 



As discussed above, the amplitude of the numerical (2, 1) 
and (4, 4) modes reaches a peak a fairly long time after 
the peak of the (2, 2) mode. Thus, in order to impose the 
first condition (in our list of five) given above, the peak of 
the EOB mode should be moved to t£ cak + Ai^ l ak . How- 
ever, the leading EOB amplitude is proportional to a 
power of the orbital frequency. This frequency decreases 
to zero at the horizon, and so the EOB amplitude drops 
to an extremely small value at t^ cak + Ai p ™ ak . By replac- 
ing v% = (Mr^fi) 2 with 1/vq, we slow the decay of these 
modes after ip eak , and can successfully move the peak of 



the mode to i£ eak 



Ai p ™ ak . This modification was also 



adopted in Ref. [28] to successfully model the (2, 1) and 
(4, 4) modes in the comparable-mass case. 



B. Comparison of calibrated EOB waveforms and 
Teukolsky waveforms for a = 

We now improve on the EOB model of Ref. [28] to more 
accurately reproduce Teukolsky waveforms. We focus on 



comparisons for the a = limit, although we discuss how 
we build our model for general spins. We start from the 
five conditions discussed at the beginning of Sec. IV A 
which allow us to compute the NQC coefficients: 



(50) 

(51) 
(52) 



d\hf°*\ 
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Tcuk/j.£m \ 
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where 



^-"pcak 



t 



tin 



t 



peak peak 
d\hj^\/dt=0 



tdQ,/dt=0 ■ 



(53) 
(54) 



(55) 



As before, the quantities hf® B given above are equivalent 
to the quantities /j^p-p 1 " 1 ^ m jr)q (10). The quantities 
on the RHSs of Eqs. (50)-(54) are now computed from 
the Teukolsky waveforms, rather than using the z/-fits of 
Ref. [28], as in Sec. IV A. Notice that Ai p ™ k in Eq. (55) 
differs from the one in Eq. (42). In fact, since the Teukol- 
sky code uses the EOB trajectory, the time difference 
between the peak of any (£, rn) mode and the peak of 
the orbital frequency is unambiguous. This was not the 
case in Ref. [28] where numerical relativity waveforms are 
used. In that case, the time difference between the peak 
of the numerical relativity (£, m) modes and the peak 
of the EOB orbital frequency depends on the alignment 
procedure between the numerical and EOB waveforms, 
and on the calibration of the EOB adjustable parame- 
ters. This is why Ref. [28] adopted the Ai p ™ k described 
in Eq. (42). 
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FIG. 10. The same as Fig. 8, but for the Teukolsky-calibrated EOB model described in Sec. IV B. 



As 



discussed 

.22 



in Sec. IV A, Ref. [28] assumed that 

t 22 



Scak = *peak and consequently set Af^ ak = 0. However, 
as seen in Tables II, III, the Teukolsky data show that 



Q 

peak 6 peak 



ij?.,, differs from zero when v = 10 3 ; this effect 



5 

is particularly pronounced for large positive Kerr spin 
parameters. The modified prescription given by Eq. (55) 
is thus quite natural. By contrast, the prescription de- 
scribed in Ref. [28] is bound to fail for all the modes in 
the test-particle limit. We point out that in Appendix 
A of Ref. [35] the authors explored how non-spinning 
waveforms can be improved, by attaching the QNMs at 
a shifted time (common to all modes) ~ 3M past the 
peak of the orbital frequency, while maintaining the de- 
termination of the NQC parameters at the peak of the 
orbital frequency. 

In the non-spinning case, we solve Eqs. (50), (51) for 



with i = 1,2,3, and set a/ m — 



= 0. In the 



spinning case, in order to not introduce spin-dependence 
at leading order, we fix a\ im and a^ lm to the values cal- 
culated for a = 0, and solve for a' ifm (with i = 3,4,5). 
As for the phase NQC coefficients, in the non-spinning 
case we solve Eqs. (52) and (53) for (with i = 1,2) 
and set b^ tm = b\ lm = 0. In the spinning case, we fix 
b\ lm and b% m to their a = values (again, in order to 
not introduce any spin-dependence at leading order) , and 
solve for b\ tm and b\ lm . 

Once the coefficients a^ em and b\ lm are known, we cal- 
culate /i™ p ~ plungc using Eq. (10) and attach the QNMs 
using Eq. (19), assuming the comb's width as in Eq. (48). 
Furthermore, we choose i^ tch = i£ oak + A^™ k with 
Ai p ™ ak given in Eq. (55). The merger-ringdown (2,2) 
mode is now attached at the time where the Teukolsky 
(2, 2) amplitude peaks, in contrast to the approach used 
in Ref. [28]. All the other merger-ringdown (£, m) modes 
are attached at the time the corresponding Teukolsky 
{I, m) amplitudes peaks, which is the same procedure 
followed in Ref. [28]. 
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■ Tcuk 
1 ,fc £m,pcak 
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■ Teuk 
^^mjpeak 


(2,2) 


-2.99 


0.001450 


-3.171xl0~ b 


0.2732 


0.005831 


(2,1) 


6.32 


0.0005199 


-7.622xl0" 7 


0.2756 


0.01096 


(3,3) 


0.52 


0.0005662 


-1.983xl0" 6 


0.4546 


0.01092 


(4,4) 


2.26 


0.0002767 


-1.213xl0" 6 


0.6347 


0.01547 



TABLE IV. Input values for the RHSs of Eqs. (50)-(54) for 
the EOB model for non-spinning black-holes used in Figs. 9 
and 10. 



Figures 9 and 10 compare these calibrated EOB mod- 
els to the Teukolsky amplitudes. Table IV lists the input 
parameters used on the RHSs of Eqs. (50)-(54). We em- 
phasize that the value of A£pg ak we reported in this table 
and that we use in our model is 2.5M smaller than the 
difference t^ cak — t^ cak . We do this because the peak of 
the Teukolsky (2, 1) amplitude is quite broad, and at the 
time ipg ak where the peak occurs the Teukolsky mode's 
frequency oscillates due to superposition of the £ = 2, 
mil modes, as discussed in Sec. HIE. (See also Ref. [34] 
and Fig. 3 of Ref. Ref. [35], which shows similar oscil- 
lations.) Although these oscillations are physical, we do 
not attempt to reproduce them in our EOB waveform, as 
their effect on the phase agreement between the EOB and 
Teukolsky waveforms is negligible. We therefore simply 
choose a slightly smaller value of Atp^ ak , ensuring that 



peak 



A/ 21 

^peak 



is within the broad peak of h^ cak . This 
in turn ensures that the oscillations in frequency do not 
impact our results. 

Figure 10 demonstrates that, by calibrating against the 
Teukolsky waveforms using the input values shown in Ta- 
ble IV, the agreement between the EOB and Teukolsky 
waveforms is considerably improved during merger and 
ringdown for the (2,1) and (4,4) modes. Improvements 
to the (2, 2) and (3, 3) modes (Fig. 9) are less significant, 
since the model of Ref. [28] works quite well for these 
modes. As we discussed in the previous section, the in- 
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FIG. 11. Comparison of Teukolsky-calibrated EOB and Teukolsky (2,2) modes for a/M = -0.5 (left panel) or a/M = 0.5 
(right panel). Upper panels show the real part of the modes; lower panels show phase and fractional amplitude differences. 



put parameters listed in Table IV are well predicted by 
the fitting formulas of Ref. [28] (see also Fig. 6). There 
is, however, noticeable improvement in phase agreement 
between Figs. 7 and 9. This is due to the fact that in the 
latter case we use the iterative procedure to compute the 
NQC coefficients aj, as discussed at the end of Sec. IV A. 

Comparing to the discussion of numerical error in 
Sec. IIIC, we note that the differences in phase and am- 
plitude between the EOB and Teukolsky modes shown 
in Figs. 9 and 10 are within the numerical errors es- 
sentially through the plunge; the differences grow larger 
than these errors during merger and ringdown. As this 
analysis was being completed, we acquired the capability 
to produce Teukolsky waveforms using the hyperboloidal 
layer method (Ref. [53]; see also Ref. [36]). We have com- 
pared phase and amplitude differences for the (2, 2) mode 
between the Teukolsky code used for the bulk of this 
analysis, and the hyperboloidal variant. We have found 
that these differences are within the errors discussed in 
Sec. IIIC. 



C. Comparisons for general spin 

We conclude our discussion of results by comparing, for 
the first time, EOB and Teukolsky coalescence waveforms 
with a^O. These waveforms are produced using the tra- 
jectory of the spinning EOB model described in Sec. II. 
As in the non-spinning case, understanding the transi- 
tion from inspiral to ringdown in the test-particle limit 
when the central black hole carries spin can help model- 
ing the plunge-merger waveforms from comparable-mass 
spinning black holes. 

As in the non-spinning case [15, 32-36, 67, 68], we ex- 
pect that the resummed-factorized energy-flux and mode 
amplitudes agree quite well with the Teukolsky data at 
least up to the ISCO, provided that the spin is not too 
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0.001542 


-1.334xl0" 6 


0.3396 


0.005095 
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0.3883 


0.004068 


0.9 


-39.09 


0.001576 


-8.102xl0~ 8 
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TABLE V. Input values for the RHSs of Eqs. (50)-(54) for 
the EOB (2, 2) mode for spinning black-holes used in Figs. 11 
(a/M = ±0.5) and 12 (a/M = 0.7 and a/M = -0.9). We 
also include data for the case a/M = 0.9, although we do not 
compare waveforms for this example. 



high. In fact, Ref. [15] showed that, in the adiabatic 
limit, the resummed-factorized (2, 2) modes agree very 
well with frequency-domain Teukolsky modes up to the 
ISCO, at least over the range -1 < a/M < 0.7. The 
relative difference between amplitudes in the two models 
is less than 0.5% when a/M < 0.5, but grows to 3.5% 
when a/M ~ 0.75. 

In this work, we focus on the (2, 2) mode comparison, 
leaving to a future publication a thorough study of the 
higher modes. We have already seen in Sec. HIE that 
as a/M — > 1, many more modes become excited during 
the plunge and merger. For this limit, the resummed- 
factorized waveforms will need to be improved in order 
to match higher-order Teukolsky modes with good preci- 
sion. 

Figures 11 and 12 compare (2,2) modes for the EOB 
and Teukolsky waves for spin values a/M = ±0.5, a/M = 
0.7 and a/M = -0.9. We build the full EOB wave- 
form following the prescription described in Sec. IV B, 
using the input parameters shown in Table V. For the 
cases a/M = 0.5 and a/M = 0.7, we also use a pseudo 
QNM (pQNM) (in addition to the standard QNMs) as 
suggested in Refs. [22, 28]. A possible physical motiva- 




FIG. 12. The same as Fig. 11, but for a/M = -0.9 (left panel) and a/M = 0.7 (right panel). 



tion of these pQNMs follows. The peak of the orbital 
frequency comes from orbits that are very close to the 
light-ring position [14], which in turn corresponds nearly 
to the peak of the effective potential for gravitational 
perturbations [69-73]. Therefore, before the orbital fre- 
quency peaks, the gravitational-wave emission is domi- 
nated by the source of the Teukolsky equations (i.e. by 
the particle); afterwards, the emission is dominated by 
the black-hole's QNMs. In the standard EOB approach, 
the waveform is a superposition of QNMs already after 
the peak of the numerical amplitude ip™ ak - However, we 
have seen that this precedes the peak of the orbital fre- 
quency by a considerable time interval: — At^ m w 12- 
40M for a/M = 0.5 and a/M = 0.7 (see Table V). To 
account for the effect of the particle emission before the 
peak of the orbital frequency, we therefore introduce a 
pQNM having frequency il>2^ NM 



2tt max (cf. Table III) 



We included this 



and decay time T22 = 
pQNM only for a/M = 0.5 and a/M = 0.7. For smaller 
spins, since Ai 22 ak is small, the pQNM would be short 
lived and would not alter our results significantly. 

As in the non-spinning case, phase and amplitude 
agreement are excellent until the ISCO. The phase dif- 
ferences remain small during the plunge, until merger, 
and grow up to ~ 0.1 rad during the ringdown. The 
amplitude difference grows to larger values, ~ 20-30% 
through merger and ringdown, because of the limitations 
of our procedure to attach the QNMs in the EOB wave- 
forms (see Ref. [28], Fig. 3 and associated discussion). 
The phase difference during the merger-ringdown for the 
case a/M = 0.7 is larger, because for larger and larger 
spins the resummed-factorized waveforms [15] perform 
less and less accurately around and beyond the ISCO. 
In the case a/M = —0.9, the disagreement between the 
EOB and Teukolsky (2, 2) becomes large and oscillatory 
during ringdown. This is a consequence of the fact that 
the oscillatory frequency behavior discussed in Sec. Ill E 
is particularly strong in this case, but we are not includ- 
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TABLE VI. The time delay At^ ak , At^ ak and At p t ak de- 
fined in Eq. (55) for a/M = -0.9, -0.5, 0.5, 0.7 and 0.9. Time 
delay information for the non-spinning case a/M = and the 
dominant (2, 2) mode are given in Table II and Table III. 



ing the associated (2, —2) QNMs in our EOB model. 

Finally, although in this paper we do not attempt to 
calibrate higher-order modes for a ^ 0, it is useful for 
ongoing work on the comparable-mass case to extract 
relevant information, such as the time delay between the 
peaks of the (£, m) modes, and the input parameters en- 
tering the RHS of Eqs. (50)-(54). In Table VI, we show 



the time delays Atp™ k ; in Fig. 13, we show 



LTcuk 

£m,peak 



and w/^peak as functions of a/M. Quadratic fits to these 



functions arc as follows: 



liTcuk 
| "22, peak | 

If, Teuk I 
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I "33, peak | 

I I, Teuk I 
| ' i 44,peak | 

Teuk 
w 22,pcak 

Teuk 
w 21,peak 

Teuk 
w 33,peak 

Teuk 
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0.001 [1.46 + 0.144 a/M + 0.00704(a/Af ) 2 ] 
0.001 [0.527 - 0.445 a/M + 0.016(a/Af) 2 ] , 
0.001 [0.566 + 0.133 a/Ml + 0.0486(a/Af) 2 ] 
0.001 [0.276 + 0.0773 a/M + 0.0405(a/Af ) 2 
0.266 + 0.129 a/M + 0.0968 {a/M) 2 , 
0.291 + 0.0454 a/M - 0.0857 (a/M) 2 , 
0.441 + 0.224 a/M + 0.163 (a/M) 2 , 
0.616 + 0.315 a/M + 0.227 (a/M) 2 . 



We postpone to future work the study 
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FIG. 13. Peak amplitude and corresponding frequency of 
Teukolsky modes as a function of spins. Top panel is the 



Tcuk 

£m,pcak I ! 



bottom is the frequency uj e 



Tcuk 



m,pcak 



peak amplitude \h 
extracted when the amplitude reaches this peak. Data is for 
the (2,2), (2,1), (3,3) and (4,4) modes, for spins a/M = 
—0.9, —0.5, 0, 0.5, 0.7 and 0.9. Solid lines are quadratic fits to 
these; the fits are given in Eq. (56). 



tTeuk 

£m,peak 



/dt 2 and dij^peak for higher-order modes. 



V. CONCLUSIONS 

The similarity of the transition from inspiral to merger 
to ringdown over all mass ratios studied in Refs. [8, 11] 
suggested the possibility of using the test-particle limit 
as a laboratory to investigate quickly and accurately the 
main features of the merger signal. The authors of Refs. 
[32, 33] were the first to exploit this possibility. They 
proposed using the EOB inspiral-plunge trajectory to 
build the source for the time-domain Regge-Wheeler- 
Zerilli equations. They also improved the EOB modeling, 
notably the energy flux and the non-quasi-circular orbit 
effects, by requiring that the EOB and RWZ leading (2, 2) 
mode agreed during plunge, merger and ringdown. 

Here, we have employed the time-domain Teukolsky 
code developed in Refs. [39 41] and extended previ- 
ous works [32, 33, 35, 36] in several directions. In 
the Schwarzschild case, we first discussed how the EOB 
model developed in Ref. [28] for comparable-mass non- 
spinning black holes performs when v = 10~ 3 for the 
leading (2, 2) mode, as well as for three sublcading 
modes, (2,1), (3,3) and (4,4). Confirming previous re- 
sults [15, 32-36, 67, 68], we found that the agreement be- 
tween the Teukolsky and EOB modes is excellent during 
the long inspiral. During the merger, whereas the agree- 
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FIG. 14. Comparison between Teukolsky-calibrated EOB 
and Teukolsky h+ and h x polarizations for a/M = 0. The 
four dominant modes (2,2), (2,1), (3,3) and (4,4) are in- 
cluded. 



ment of the (2, 2) and (3, 3) modes is still good, that of the 
(4, 4) and (2, 1) is not very satisfactory. We find that this 
is due to the irregular behavior of the numerical-relativity 
input values for the peak of the mode amplitude and the 
gravitational frequency at that peak. This motivates the 
need for more accurate numerical-relativity data for these 
higher-order modes, which will presumably be available 
in the future. By calibrating the EOB model using input 
values directly extracted from the Teukolsky modes (Ta- 
bles II and III), we found very good agreement for the 
four largest modes. In Fig. 14, we compare h + and h x 
constructed for these four modes, using 

h+(0, <(>, t) - ih x (9, cf>,t) = J2 -2Y em (0, 4>) h lm {t) . (57) 

l,m 

The sum here is over (£,m) = (2, ±2), (2, ±1), (3, ±3) 
and (4, ±4). The agreement between EOB and Teukolsky 
polarizations is very good as expected. There are some 
minor differences during the ringdown, which are mainly 
due to the underestimated ringdown amplitudes of the 
(2, 2) and (3, 3) modes in the EOB model. 

Moreover, for the first time, we employed the EOB 
inspiral-plunge trajectory to produce merger waveforms 
for quasi-circular, equatorial inspiral in the Kerr space- 
time. The energy flux in the EOB equations of motion 
uses the factorized resummed waveforms of Refs. [13, 15]. 
We calibrated the leading EOB (2, 2) mode for spins 
a/M = —0.9, —0.5, 0.5, 0.7, and extracted information on 
the subleading modes. We also investigated the high spin 
case a/M = 0.9. We found that several modes which are 
subleading during the inspiral become relevant during 
plunge and merger. The major new feature of the EOB 
calibration (based on Teukolsky data) is that we relaxed 
the assumption used in previous papers [18-20, 22-28] 
that the matching of the QNMs for the leading (2, 2) 
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mode occurs at the peak of the orbital frequency. In 
fact, we found that the peak of the orbital frequency does 
not occur at the same time as the peak of the Teukolsky 
(2, 2) mode, and that the time difference grows as the 
spin parameter increases. Our work represents a first 
step in exploring and taking advantage of test-particle 
limit results to build a better spin EOB model in the 
comparable-mass case [45]. 

In the future, we plan to extend this work in at least 
two directions. First, we want to calibrate the EOB 
model in the test-particle limit for higher spins and for 
higher-order modes, and to connect it to the spin EOB 
model in the comparable-mass case [23, 45]. To achieve 
this goal, we would need to introduce adjustable param- 
eters in the functions pi m in Eq. (11) to improve the 
resummed-factorized energy flux and amplitude modes 
for large spin values. 

In our future analyses, we will use a Teukolsky code 
which uses hyperboloidal slicing [36, 53]. Although we 
were able to achieve similar accuracy by extrapolating 
our results from finite radius to future null infinity, hy- 
perboloidal slicing is far faster, and has proven to be very 
robust. Second, we would like to extend this model to 
inclined orbits. To tackle this case, we need to generalize 
the resummed-factorized waveforms to generic spin orien- 
tations. If we were only interested in extracting the input 
values, as in Tables IV, V, it might be sufficient to use 
the hybrid method suggested in Ref. [48]. In this case, 



we could use in the EOB equations of motion the energy 
flux computed with a frequency-domain Teukolsky code 
[61], but extend it to plunging trajectories. 

Finally, besides improving the EOB model, the possi- 
bility of generating quickly and accurately merger wave- 
forms in the test-particle limit will allow us to investi- 
gate several interesting phenomena, such as the distribu- 
tion of kick velocities for spinning black- hole mergers [41], 
the energy and angular-momentum released when a test 
particle plunges into a Kerr black hole [74, 75], and the 
generic ringdown frequencies suggested in Refs. [49, 50]. 
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